function tests = testLPF()% {{{
    close all;
    tests = functiontests(localfunctions());
end % function }}}
function setup(testCase)% {{{
    % 设计一个低通滤波器
    fs = 1000;  % 采样频率1kHz
    T = 1;      % 持续时间 1s 
    dt = 1/fs;
    N = round(fs * T);
    t = 0:dt:T-dt; % 时间序列
    f1 = 5;    % 50Hz
    s1 = sin(2*pi*f1*t); % 50Hz正弦波信号
    a1 = 0;
    f2 = 40;   % 200Hz
    s2 = sin(2*pi*f2*t); % 200Hz正弦波噪声
    sMix = a1*s1 + s2; % 合成信号

    n = 50;     % 设计30阶低通滤波器
    wCutOff = 0.25; 
    fCutOff = 250;
    b = fir1(n, fCutOff/fs); 
    testCase.TestData = struct('Time', t, 'InputSignal', sMix, 'FIRFilterCoeff', b, 'FIRFilterOrder', n, 'InputFrequencies', [f1, f2], 'SampleFrequency', fs);
end % function }}}
function testFilter(testCase)% {{{
    t = testCase.TestData.Time;
    sMix = testCase.TestData.InputSignal;
    b = testCase.TestData.FIRFilterCoeff;
    y1 = filter(b, 1, sMix); % 应用滤波器
    % 绘制结果
    figure('WindowStyle', 'docked');
    subplot(2, 1, 1);
    plot(t, sMix);
    title('原始信号');

    subplot(2, 1, 2);
    plot(t, y1);
    title('滤波后信号');
end % function
function testFreqz(testCase)
    b = testCase.TestData.FIRFilterCoeff;
    n = testCase.TestData.FIRFilterOrder;
    freq = testCase.TestData.InputFrequencies;
    t = testCase.TestData.Time;
    fs = testCase.TestData.SampleFrequency;
    f1 = freq(1);
    f2 = freq(2);
    M = numel(t);
    inputSignal = randn(1, M) + 1j*randn(1, M);
    [inputSpectrum, frequency1] = fourierTransform(inputSignal, t);
    outputSignal = filter(b, 1, inputSignal);
    [outputSpectrum, frequency2] = fourierTransform(outputSignal, t);
    assert(all(frequency1 == frequency2));
    [H, f] = freqz(b, 1, frequency1, fs);
    w1 = (0:M-1)'/M*pi;
    H1 = abs(sum(b .* exp(1j*w1.*(0:n)), 2));
    w = 2*pi*f/fs;
    % assert(max(abs(w - w1)) < 1e-15);
    % assert(max(abs(H) - abs(H1)) < 2e-15);
    % plot(w, abs(H), '.', w, abs(H1), 'o');
    % Hmin = min(db(H));
    % Hmax = max(db(H));
    outputSpectrum1 = inputSpectrum .* conj(H);
    max(abs(abs(outputSpectrum1) - abs(outputSpectrum)))
    figure('WindowStyle', 'docked');
    plot(f, db(H), frequency1, db(inputSpectrum), frequency2, db(outputSpectrum));
end % function}}}

function testFilterCheck(testCase)
    b = testCase.TestData.FIRFilterCoeff;
    y1 = testCase.TestData.InputSignal;
    t = testCase.TestData.Time;
    n = testCase.TestData.FIRFilterOrder;
    % test filter
    y2 = filter(b, 1, y1);
    y3 = zeros(size(y2));
    for k = 1:numel(y3)
        y3(k) = 0;
        for i = 0:min(k-1,n)
            y3(k) = y3(k) + b(i+1) * y1(k-i);
        end % for
        % assert(abs(y3(k) - y2(k)) < 1e-15);
    end % for
    figure('WindowStyle', 'docked');
    plot(t, y1, 'b-', t, y2, 'r.-');
end % function}}}

